﻿class verlet
{
private:
	caja *box;
	///Paso de tiempo Dt
	double Dt;
	///Dt2
	double Dt2;
	///Dt/2
	double Dth;
	///Dt2/0.5
	double Dt2h;
	///Uso del termostato
	bool termo;
	///Coeficiente de friccion
	double xgamma;
	///sigma = 2*gamma*mi*kB*T
	double sigma;
	///C1 = Dt*(1-1/2*Dt*gamma)
	double C1;
	///C3=0.5*sigma*Dt
	double C3;
	///D1=1/2*Dt*(1-gamma*Dt)
	double D1;
	///D3=gamma*C1*mi
	double D3;
	///D4=sigma*(1-1/2*gamma*Dt)
	double D4;
	///Para el wiener process Wien=sqrt(3*Dt)
	double wien;
	///Para el motor Verlet-Langevin
	double temp, diff;
    double xs, xsd;
    double c0,c1,c2,c1dt,c2dt,c2dt2,c1mc2dt;
    double sigmaR, sigmaV;
    double corrrV,umcorrV2SQ;
    double **drg, **dvg;
    bool langevin;

public:

	verlet(caja *bx, double xDt, bool xtermo)
	{
		box=bx;
		Dt=xDt;
		Dt2=Dt*Dt;
		Dth=Dt*0.5;
		Dt2h=0.5*Dt*Dt;
		termo=xtermo;
		langevin=false;
	}

    verlet(caja *bx, double xDt, double xtemp, double xdiff)
    {
        box=bx;
        temp=xtemp;
        diff=xdiff;
        Dt=xDt;
        xs=kB*temp/diff;
        xsd=xs*Dt;  //(kBT·t/(m·Dif))
        c0=exp(-xsd);
        c1=(1.0-c0)/xsd;
        c2=(1.0-c1)/xsd;
        c1dt=c1*Dt;
        c2dt=c2*Dt;
        c2dt2=c2*Dt*Dt;
        c1mc2dt=(c1-c2)*Dt;
        sigmaR= sqrt(Dt*Dt*kB*temp/xsd*(2.0-(3.0-4.0*exp(-xsd)+exp(-2.0*xsd))/xsd));
        sigmaV= sqrt(kB*temp*(1.0-exp(-2.0*xsd)));
        corrrV=Dt*kB*temp/xsd*(1.0-exp(-xsd))*(1.0-exp(-xsd))/(sigmaR*sigmaV);
        umcorrV2SQ=sqrt(1.0-corrrV*corrrV);

        drg = new double*[box->npart];
        dvg= new double *[box->npart];
        for (int i=0;i<box->npart;i++)
        {
            drg[i]=new double[3];
            dvg[i]=new double[3];
        }
       // dtem  = dtem*ldtem*sign(eins,etemp-atemp);
    }

	verlet(caja *bx, double xDt, double sgamma, bool xtermo)
	{
	    langevin=false;
		xgamma=10.0;
		box=bx;
		Dt=xDt;
		Dt2=Dt*Dt;
		Dth=Dt*0.5;
		Dt2h=0.5*Dt*Dt;
		termo=xtermo;
		sigma = sqrt(2.0*xgamma*box->temp0);
		C1 = Dt*(1.0-0.5*Dt*xgamma);
		C3=0.5*sigma*Dt;
		D1=0.5*Dt*(1.0-xgamma*Dt);
		D3=xgamma*C1;
		D4=sigma*(1.0-0.5*xgamma*Dt);
		wien=sqrt(3.0*Dt);

	}

	///Motor LeapFrog con link cell list
	void  LeapFrogLinkCellList(bool calcSv);
	//Motor Verlet clasico con link cell list
	void  verletLinkCellList(bool calcSv);
	///Motor dinamica browniana de Paul y Yoon (PRE 52,2,2076 1995)
	inline void PaulYoonEngine();
	///Wiener process
	double wiener(void);
	///Copia literal del motor de Juan Pablo
	void posvel(void);
	///Motor dinamica browniana de Paul y Yoon (PRE 52,2,2076 1995) con linkCell list
	void PaulYoonEngineLinkCell(bool calcSv);
	void LeapFrogSimple(bool calcSv);
	///Inicializacion del motor Verlet-Langevin
	int initVerletLangevin(void);
	///Motor Verlet-Langevin
	void  verletLangevin(void);
	///Generador variables aleatorias gaussianas
	double gaussianRand(void);
	///Componentes aleatorias posicion-velocidad
	void gaussian(void);


};

void  verlet::LeapFrogSimple(bool calcSv)
{
	///si calcSv = true, se calcula el 'shear stress tensor'
	double dv[3];
	double parke=0.;
	///Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz;

	///Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for(int i=0;i<box->npart;i++)
	{
		for (int j=0;j<3;j++)
		{
			box->part[i].pos[j]+=(box->part[i].vel[j])*Dt+(box->part[i].force[j])*Dt2h;
			///Calculamos la velocidad en t+Dt/2
			if(termo)
			{
				///Empleamos un reescalado de las velocidades como termostato
				box->part[i].vel[j]=(box->part[i].vel[j])*sqrt(box->temp0/box->temp)+Dth*box->part[i].force[j];
			}
			else
			{
				///No utilizamos termostato
				box->part[i].vel[j]+=Dth*box->part[i].force[j];
			}
		}
		box->boundary(i);
	}
	///Calculamos las fuerzas, epotencial y el virial
	interaccionSimple(box);

	///Velocidad en t+Dt
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]+=(box->part[i].force[j])*Dth;
			///Actualizamos la componenete de la energia cinetica
			//parke+=SQR(box->part[i].vel[j]);
			///Calculamos la deriva media de las velocidades de la caja
			dv[j]+=box->part[i].vel[j];
		}

	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	///Eliminamos la deriva de las velocidades de las particulas
	for(int i=0;i<box->npart;i++)
	{
		for(int j=0;j<3;j++)
		{
			box->part[i].vel[j]-=box->vel[j];
			///Actualizamos la componenete de la energia cinetica
			parke+=SQR(box->part[i].vel[j]);
		}
		if(calcSv)
		{
			///Calculamos la componente de velocidad del tensor de esfuerzos
			svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
			svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
			svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
		}
	}

	if(calcSv)
	{
		///Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy/box->vol;
		box->svxz=svxz/box->vol;
		box->svyz=svyz/box->vol;
	}

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	//Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol; //IMPORTANTE-> Falta introducir la correccion del cutoff a la presion
}

inline void verlet::PaulYoonEngine()
{
	///Motor de Paul con calculo de interacciones entre todas las particulas de la caja
	///(sin cutoff). Comprobar para distintos numeros de particulas
	double mass=1.0;

	///Velocidad media de las particulas de la caja
	double dv[3];
	double xwien;
	///Elementos del 'microscopic stress tensor' asociado a la componente cinetica
	double svxy, svxz, svyz, svzz;
	///Componente de la energia cinetica
	double parke=0.;
	///Inicializamos los elementos del tensor de esfuerzos
	svxy=svxz=svyz=svzz=0.0;
	///Para calcular la velocidad media del sistema
	dv[0]=dv[1]=dv[2]=0.0;

	for (int i=0;i<box->npart;i++) 	{
		for(int j=0;j<3;j++){
			xwien = wiener();
			double pp = box->part[i].vel[j];
			box->part[i].pos[j]+=C1*(box->part[i].vel[j])+Dt2h*box->part[i].force[j]+C3*xwien;
			xwien = wiener();
			box->part[i].vel[j]+=D1*(box->part[i].force[j])-D3*(box->part[i].vel[j])+D4*xwien/mass;
		}

		#ifdef FLUJOS
		///Pruebas con flujos
		    if(box->part[i].pos[2]+0.5*box->lbox[2]<1.){
                box->part[i].vel[1]=0.;
		    } else {
                box->part[i].vel[1]+=box->vel_shear*(box->part[i].pos[2]/box->lbox[2]+0.5);
		    }
		#endif
		///Condiciones de contorno periodicas
		box->boundary(i);
	}
	///Calculamos las fuerzas en la nueva posicion
	#ifdef CONFINED
	interaccionConfined(box);
	#endif
	#ifndef CONFINED
	interaccionSimple(box);
	#endif

	///Acualizamos la componente de las velocidades que depende de la fuerza en el nuevo paso de tiempo
	for (int i=0;i<box->npart;i++) 	{
		for(int j=0;j<3;j++){
			box->part[i].vel[j]+=Dth*(box->part[i].force[j]);
			dv[j]+=box->part[i].vel[j];
		}
	}

	///Calculamos la velocidad media de la caja
	for(int i=0;i<3;i++)
		{box->vel[i]=dv[i]/double(box->npart);}

	///Eliminamos la velocidad de deriva del sistema
	for (int i=0;i<box->npart;i++){
		for(int j=0;j<3;j++) {
			box->part[i].vel[j]-=box->vel[j];
			parke+=box->part[i].mass*SQR(box->part[i].vel[j]);
		}
		#ifdef CALC_STRESS
			///Calculamos la componente de velocidad del tensor de esfuerzos
			#ifndef POLYDISPERSE
                svxy+=(box->part[i].vel[0])*(box->part[i].vel[1]);
                svxz+=(box->part[i].vel[0])*(box->part[i].vel[2]);
                svyz+=(box->part[i].vel[1])*(box->part[i].vel[2]);
                svzz+=(box->part[i].vel[2])*(box->part[i].vel[2]);
			#endif

			#ifdef POLYDISPERSE
                svxy+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[1]);
                svxz+=(box->part[i].mass)*(box->part[i].vel[0])*(box->part[i].vel[2]);
                svyz+=(box->part[i].mass)*(box->part[i].vel[1])*(box->part[i].vel[2]);
                svzz+=(box->part[i].mass)*(box->part[i].vel[2])*(box->part[i].vel[2]);
			#endif
		#endif
	}

	#ifdef CALC_STRESS
		///Actualizamos el valor de la componente de velocidad del tensor de esfuerzos
		box->svxy=svxy/box->vol;
		box->svxz=svxz/box->vol;
		box->svyz=svyz/box->vol;
		box->svzz=svzz/box->vol;
	#endif

	///Introducimos la energia cinetica del sistema
	box->ekin=0.5*parke/double(box->npart);
	///Calculamos la temperatura de la caja
	box->temp=parke/3.0/double(box->npart);
	///Calculamos la presion de la caja
	box->presion=box->rho*box->temp+box->virial/box->vol;
}


double verlet::wiener(void)
{
	double alea = double(rand())/RAND_MAX;
	double dd=0.02236068;
	return 1.732050808*dd*int(3.*alea-1.5);
}
